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SUMMARY 


Error analyses of di stant-satel 1 i te- to-close-satel 1 i te 
range-rate, satel 1 i te-to-sea altimetry, and ground station to 
satellite range are made by simulations in which observational 
variances are assumed, observation equations are formed, and 
normal equations incremented. The final normal equation matrix 
is inverted to obtain standard deviations and correlation coefficients. 

The natural parameters solved for are the broad variations 
of the gravity field, represented by harmonic coefficients; 
local variations of gravity, represented by point masses; and 
the departure of the sea level from the geoid, represented by 
area means. 

A standard case of a low (263 km) polar close satellite, 
three equatorial geosynchronous satellites, and eight ground 
tracking stations is set up. Sigmas of +0.5 mm for range-rate, 

+0.5 m for range, and +2 m for altimetry are assumed. Twenty- 
one variations from this standard case are tested. 

The principal conclusions: 

1. Given +0.5 mm/sec range-rate tracking from distant 
satellites and +2 m altimetry, the range-rate is considerably 
more effective in determining variations of the gravity field. 

2. For a satellite altitude of 260 km and +0.5 mm/sec 
range-rate, the average resolution of determination of gravity 
variations will be about 3°, or 350 km. 

3. The accuracy is sensitive to satellite altitude: a 

change of 60 km makes a difference of a factor of 50% in accu- 
racy of determination of a given size element in the gravity 
field. 
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4. The accuracy is sensitive to the number of geosynchronous 
satellites from which the ciose satellite is tracked; apparently 
it is advantageous for the effect of a given element in the 
gravity field to be sensed from more than one direction. 

5. A second close satellite at a lower inclination (say 
45°) in addition to a polar satellite improves the accuracy 
by a factor of about 30%. 

6. The accuracy of determination of the gravity varia- 
tions is insensitive to the number and distribution of ground 
tracking stations, the accuracy of tracking station location 

and even the accuracy of the distant satellite orbit determination. 

\ 

7. Altimeter accuracy of +20 cm or better is needed to 


determine departures of mean sea level from the geoid. 
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I INTRODUCTION 


Purpose 

The purpose of this work is to estimate the accuracies with 
which earth parameters can be determined by new tracking tech- 
niques for given tracking station and satellite orbit configura- 
tions. It therefore is intended to aid in the optimization of 
satellite orbit specifications, tracking station locations, etc. 

The tracking techniques considered in this report are: 
satel 1 i te- to-sate 1 1 i te range-rate; satel 1 i te- to-sea altimetry; 
and ground stat i on- to-satel 1 i te range. There is no limitation 
on the satellite orbits considered, other than that, in a satellite 
to-satel lite tracking pair, one satellite is assumed to be sig- 
nificantly perturbed, while the other is not: i.e., one is much 

lower than the other. 


Method 

The method employed is to assume standard deviations for 
the tracking, construct a hypothetical series of observations, 
form observation equations and thence normal equations, and in- 
vert the final normal equation coefficients to obtain a covariance 
matrix for the parameters. Since the observations are hypothetical 
rather than real, for purposes of the error analysis orbits can 
be assumed to be secularly precessing ellipses, and the pertur- 
bations used to determine earth parameters linearly superimposed 
thereon. 
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Allowance is made in the procedure for a priori sigmas (stan- 
dard deviations) for the parameters. These sigmas are estimates 
of the present accuracies to which the parameters are known. 

If the systems analyzed can make a significant contribution to 
the determination of a parameter, then the sigma obtained by 
taking the square root of the final covariance matrix should 
be appreciably less than the a priori sigma. While this pro- 
cedure may seem less ideal than assuming no a priori sigmas what- 
soever, it has the practical advantage of avoiding computational 
difficulties with ill-conditioned matrices. 

Provision was made for the modeling of biases and other 
systematic errors. However, experience with such modeling to 
an elaborate extent in error analyses of lunar ranging under 
grant NGL 05 - 007-283 [ Kaula . 1972 ] indicated that the results 
still give optimistically low sigmas. Hence such an effort 
is not made in the present analysis. The absolute sigmas in 
the results are therefore not to be accepted at face value. 

The changes in sigmas with changes in configuration are prob- 
ably quite meaningful. Hence, given tfh-at the tracking results 
in significant improvement from rather simple considerations, 
the error analysis aids in choosing the optimum configuration. 

I I SYSTEM COMPONENTS 


Locat i ons 


The system which the computer programs can analyze includes 



four types of sites for tracking instrumentation: 


1) ground tracking stations, 

2) close artificial satellites, 

3) distant artificial satellites, 

4) retro- ref 1 ectors on the moon. 

The functional distinction between a "close" and a "distant" 
artificial satellite is that the former is measurably perturbed 
by longitudinal variations in the earth's gravitational field, 
while the latter is not. 

Parameters pertaining to a ground station location are the 
three components of its position. Parameters pertaining to the 
satellites are their constants-of- i ntegrat i on at epoch. 

T rack ? nq 

Between the four types of locations listed above, three 
types of tracking were assumed: 

1) range from ground tracking stations to: 

close satel 1 i tes , 

d i stant satel 1 i tes , 

retro- ref lectors on the moon; 

2) range-rate from distant satellites to close satellites; 

3) altimetry from close satellites to the sea surface. 
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The accuracies usually assumed for the tracking types were 
those felt to be feasible with present technology at the 1969 
Williamstown study: +0.5 m for ground- to- satel 1 i te range, +0.5 mm/sec 

for satel 1 i te-to-satel 1 i te range-rate, and +2 m for radar alti- 
metry. For one run, accuracies of +0.05 mm/sec and +0.10 m were 
assumed for the latter two types [ MIT . 1970], 

Although ranging to retro- ref lectors on the moon is included 
in the capability of the program, the analyses described below 
do not include it because of the emphasis on determination of 
the gravi tational field and mean sea level. 

Consideration was also given to the inclusion of radio 
interferometry (VLB! ) among the tracking types. An interfero- 
meter measures the difference in time of receipt of a signal 
at two stations. Hence, for objects at distances not extra- 
ordinarily larger than the baseline length the observation is 
best represented as a difference in the ranges from the two 
baseline ends. As such, the VLBI is the same as coupled range 
devices, with consequently fewer tracking opportunities. Hence 
for a given accuracy the error analysis will show the VLBI to 

i 

be poorer than ordinary ranging. 

Whether the VLBI can attain higher accuracy depends, of 
course, on atmospheric refractional effects. These will be 
larger than for laser ranging because of water vapor and higher 
zenith angles. It was concluded that a meaningful comparison 
of VLBI and ranging in an error analysis would first require 
a detailed study of the relative atmospheric refraction effects, 
hence VLBI was not included in the studies described in this 
report. 
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I I I PARAMETERS 


Class I f icat ion 

The natural parameters the program is designed to test 
the system for ability to determine can be grouped into four 
classes: 

1) Earth-geometr ical : rotation, polar wobble, station 

position and drift. 

2) Earth-long periodic gravitational: tides, mass shifts, 

zonal harmonics. 

3) Earth-short periodic gravitational: spherical harmonic 

coefficients, mass distributions, mean sea level. 

4) Moon -geometr ical : lunar ephemaris, retro- ref 1 ector 

locations, physical librations. 

Analyses pertaining to class k) were carried out under grant 
NGL 05-007-283 (Error Analysis of Lunar Ranging) [ Kaula . 1972]. 

Some work was done on class 2) in furtherance of the objectives 
of grant NSR 005-020-379 (Definition of a Drag-Free Satellite 
for Geodynamics) [ Stanford . 1970]. An error analysis was made 
in which the parameters determined from satellite perturbations 
were the fixed zonal harmonics of the gravi tational field, north- 
south variations in the tidal Love numbers and phase lags, and 
north-south seasonal mass shifts. The results were somewhat 
discouraging; for variations of wavelength shorter than representable 



7 


by a 2nd degree harmonic, the sigmas were larger than the physi 
cally plausible magnitudes. 

Definition and Form 


The emphasis of this study is therefore on the determination 
of parameters in class 3), short periodic gravi tat ional , with 
some investigation of their interaction with station position, 
in class 1). The functional distinction between "spherical har- 
monic coefficients" and "mass distributions" in class 3) is that 
the former represent variations in the gravitational field of 
wavelengths long enough that their integrated effect on the satel- 
lite position is large enough to be determined by tracking from 
ground stations alone, while the latter require either satellite- 
to-satellite range- rate or altimetry to be determinable. Mean 
sea level is here taken to be the displacement of the sea level 
from the geoid due to temperature, salinity, and steady-state 
air pressure and wind drag. 

To economize in computer storage, the spherical harmonics 
considered were limited to the highest degrees which are well- 
determined by "classical" satellite geodesy techniques f Gaposchk i n 
& Lambeck . 1971]. These we estimate to be the 11th and 12th 
degrees. Furthermore, only a limited number of orders n were 
taken. The total "long wavelength" gravity field was thus as- 
sumed to be represented by harmonics l,m = 11,0; 11,3; 11,6; 

11,9; 12,0; 12,3; 12,6; 12,9; and 12,12. 

To economize in partial derivative calculation, the shorter 
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wavelength mass distribution was assumed to be represented by 
point masses at uniform intervals. These intervals are neces- 
sarily less than 180°/12*15° if the point masses are to be non- 
zero. To economize in computer storage, these point masses were 
assumed to exist only within a test area covering a limited por- 
tion of the earth's surface: 30° x 30° , so as to encompass two 

wavelengths of the spherical harmonics. To further economize, 
on computer runs testing aspects other than how finely the gravity 
field can be resolved, the spacing was assumed to be 10°. 

A Priori Sigmas 

On those runs where interaction with tracking station position 
was taken into account, the error in each coordinate of the sta- 
tion position was assumed to be +30 m, probably triple the actual 
uncertainty. 

The a priori sigmas for the spherical harmonic coefficients 
were taken to be those given by the 10 rule-of-thumb. 

The a priori sigma for a mass point was obtained as fol- 
lows. Given the rms coefficient of potential, 6^(V), the power 
spectrum of the potential is 

0*(V) = (24+1 (1) 
and of a surface mass distribution p 


o?(p) * (24+1 ) 3 M4) 2 «*(v) 


( 2 ) 
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[ Kaula . 1968, p. 67]. 

2 

The mean square total anomalous mass a s (m) for a square of side 
1 ength swill be 

tt/s 

Og(m) = s** ^ a£(p) (3) 

1=2 


The amount residual to a specific harmonic degree t x therefore 
will be 


A a s (m) 


rr/s 


I 


'-"V 1 


(4) 


Numerical values for various size squares, using l » 12 

/x 


s 

Ao s (m) 

2.5° 

0. 15 x 10" 8 

5° 

0.60 x 10' 8 

10° 

0.5 x IQ’ 7 

15° 

0.0 


The unit of ACT g (m) is the earth's mass. 

The a priori sigma used for the departure of mean sea level 
from the geo id is based on calculation of the mean sea level 
from temperature and salinity data l Stommel . 1966, p. 180], 

This appears to be about +20 cm. This value was used in the 
analyses; it is probably an underestimate, but the results obtained 
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did not Indicate anything other than the obvious conclusion that 
an oceanographically valuable altimeter should have an accuracy 
small compared to the expected variations. 

Partial Derivatives 

The Kepler elements of the orbit were used as an intermediary 
in writing the partial derivatives with respect to spherical 
harmonic coefficients and the orbital constants-of- i ntegrat ion . 

The procedures in Chapter 4 of Kaula [1966] were followed for 
all calculations relating Kepler elements to Cartesian coordinates, 
ranges, etc. The principal new development related to the range- 
rates. To keep the partial derivatives with respect to the mass 
points simple, the satel 1 i te- to-satel 1 i te Doppler was treated 
as an acceleration : i.e., range-rates r at time intervals fit 

are replaced by the acceleration 

*r = r/6 1 (5) 

For a mass point spacing s (angular), the time interval fit is 
taken to be s/n, where n is the mean motion of the satellite. 

Given distant satellite inertial coordinates X^, close satel- 
lite inertial coordinates * , and mass point body fixed coordinates 

u , the mass point coordinates must first be converted to in- 
~m 

ertial by applying the negative (clockwise) rotation through 
the Greenwich Sidereal Time 0 



( 6 ) 




using the convention that the rotation axis is the 3-axis. Then 
letting r^, r cm , r dm be the three distances calculated from the 
inertial coordinates, the partial derivative will be for the com 
ponent of an acceleration along the line cm in the direction dc 



using the cosine law to calculate the angle at the close satel- 
lite, ty c . The negative sign appears because a positive accelera- 
tion is counted to be away from the distant satellite. Alter- 
natively, cos$ can be expressed in terms of the vector dot product 


• • 
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( 8 ) 


In (8) and hereafter, summation over subscripts repeated in a 
product, i, j, or k, applies. 

• • 

To obtain the component of acceleration r^ along the line- 
of- sight between the two satellites, define 
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Then 


U\ = *ci " *di 


r i± 

r i " 


(9) 


( 10 ) 


» •• v * * • / 

\t = U-tiX-ti + *iiUi)J r i " r iniH\/ r 




( 11 ) 


2/ r 3 


( 12 ) 


For the partial derivatives of the acceleration (12) with 
respect to the constants-of- i ntegrat ion and the spherical har- 
monic coefficients, partial derivatives with respect to oscu- 
lating Kepler elements s c j and s^j are needed 


• • 






ax 

"as 


c i 
cj 


( 13 ) 
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where, from (10) 


dr 


l X/i . sx ci 


as cj 


ss cj 


(14) 


* 

The partial derivatives dX . /as . and dX ./ ds ; are given 

C I I CJ C 1/ Cl 

by Kaula [ 1966, pp, 67-68] 


dX ci 


dR : 


iJ 


d(o, T, to ) d(n, i, a) ) q j 


dx ci 


d(a, e, M) R ij d(a, e, M) 


dq ; 


(15) 


**ci 


>R. 


d(n, i , «)) d(n, i , w) q j 


dx • 

*ci 


d (a , e, M) R i j d(a, e, M) 


aq j 


(16) 


where Rjj is an element of the rotation matrix R Xq [ Kaula . 1966, 
P- 18] 


-Xq 


R q = R 3 ( - Ci ) R , (- I ) R 3 (-u)) 


(17) 


and qj , qj are components of position and velocity in a coordinate 
system with the 3-axis normal to the osculating orbital plane 
and the 1-axis toward perigee [ Kaula . 1966, p, 24] 


/* a (cos E - e ) 

a a f- 


e sin E 
0 


(18) 
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Simi larly, 


( 


sin E 


r, 2 n V na 

'1-e cos En - r 

c ‘" i - e cos E 


0 


(19) 


• * 

ax. 


aR, 


'ci 0 I iJ * # 
a (o, i , oj ) a (n, I , tu)"j 

^Ci _ n 8q j 

a(a, e, M) i j a(a, e, M) 


( 20 ) 


where the acceleration vector q can be taken as derived from 

2 


the central term -\sjr‘ 


9. ■ -uaA 3 


Bq 

as 




aq 

*s. 


m 


( 21 ) 


The notation used here for Kepler elements is the same as in 
Kaul a f 19661 : Q, I , w are, respectively, nodal longitude, in- 

clination, and argument of perigee; a, e, M are semi-major axis, 
eccentricity, and «mean anomaly; n, E are mean motion and eccen- 
tric anomaly. The partial derivatives of the Kepler elements 
with respect to the orbital constants-of- i ntegrat ion and spheri- 
cal harmonic coefficients of the gravity field are given in Kaul a 
[1966, pp. 69-70]. 

The other new development not discussed in Chapter 4 of 
Kaula [ 19661 is radar altimetry. For the satellite altitude 
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above the sea we can write 

h = r - (r e (<p) + N + v) (22 ) 

where r is the orbital radius, r & is the radius of the refer- 
ence ellipsoid, N is the geoid height, and v is the departure 
of the sea level from the geoid. For the purposes of an error 
analysis we can assume the flattening of the reference ellip- 
soid to be negligible, and thus concentrate on the quantities 
dependent on significantly unknown parameters. These include 
the orbital constants-of- i ntegrat i on , s^; the long wavelength 
variations of the gravitational field, expressed by the harmonic 
ceofficients C^ m , S^ m ; the short wavelength variations of the 
gravi tational field, expressed by the mass points itu ; the de- 
parture of the tide height from the equilibrium value, 6?; and 
the departure ri of sea level from the tide height because of tern 
perature, salinity, wind stress, and short term pressure varia- 
tions. Thence 


r = 



dr 


as 


■6 S 


Ol 


Ol 



+ 6r(m) 


(23) 


assuming that the direct and equilibrium tide effects on the 
orbit are adequately known, and that the tidal 6? and oceano- 
graphic perturbations r\ are of negligible effect. These assumptions 
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are all quite safe, at least for orbit durations of a few days. 

A larger effect will be 6r(m), the perturbation of the orbital 
radius by the mass points representing the shorter wavelength 
gravi tat ional variations. However, it is awkward to calculate 
6r(m), and furthermore it is certainly smaller than the mass 
point effect on the geoid height, 6N(m), or on the range-rate 
6r^(m). Hence we assume it also to be negligible. For the 
partial derivatives of the radius with respect to the constants- 
of- integration and spherical harmonic coefficients, we again 
use the osculating Kepler elements as intermediaries. From 

r = a{1 - e cos E) (24) 

we get 

a (a, e, K) * (a* sln E " a cos E ' ^ rC sin ( 2 5) 


The geoid height N at cpj,Xj can be written as 


N = 


R I P tm( sin ”l)[ C tm cos "*| +S im sin mX ij 

K* i ni 


9^ r ij 


(26) 


where r.j is the distance fromcp.,Xj to the mass point location, ’ 
R is the radius of the earth, and g is the gravitational ac- 
celeration. Thence 
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RP lmi[ cos m *i * 


sin mX 


*1 

iJ 





( 27 ) 

( 28 ) 


A final assumption made was that the number of passes over 
each square of sea surface was sufficient to discriminate its 
mean departure from an equipotential , ari 0 , from the time varying 
components of v. 

To summarize, the partial derivatives pertaining to the 
altitude over a subsatellite point location cp.,Xj are 


as oi 


as 


95 i as oj 


_£h 


JUL. 
as 


a ( C -tm' S ^m) SS ' 3 ( C ^,m' S tm) 


as . _ r 

“ * ~TZ — r - R p ,J cos mx. , sin mX ; I 

i af F ' ^ 1 


( 29 ) 


( 30 ) 


J*L R_ j x : 

am. Mr . ’ ^ J 


( 31 ) 


J 


>J 


-2S1. I * j 


( 32 ) 


where s is point spacing. 


—an. = 

aT >oj 


- 6 j j , the Kronecker delta, 


(33) 


Standard Case and Variations 


To keep computer time within reasonable bounds, and to 



18 


minimize confusion, it is desirable (as in most error analyses) 
to choose one set of parameter values to define a standard con- 
figuration. The effects of changing parameters are then tested 
by one-at-a- time variations from the standard case. 

The parameters can be classified into four groups for gene- 
ral discussion. We give here only the significant parameters 
and their values. The full set required for a computer run 
are best learned from the source program and sample input avail- 
able on request. 

1 . Stations 


The number, capabilities, percent clear weather, and 
locations of ground tracking stations can be specified. For 
the standard case, we assumed eight wel 1 -di str ibuted stations 
with capability of laser ranging to both close and distant satel 
1 i tes and 50% probability (random day-to-day) of clear weather. 


The locations: 


38° N, 23° E: 

Greece 

35° N, 139° E: 

Japan 

20° N, 204° E: 

Hawai i 

35° N, 282° E: 

North Carolina 

33° S, 18° E: 

Cape of Good Hope 

3 1° s , 136° E: 

South Australia 

26° S , 2 16° E : 

Rapa (south of Tahiti) 

33° S, 290°E: 

Central Chile 
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2. S ate 1 1 i tes 

The number and orbital elements of each type of satel- 
1 i te can be speci f ied. For a close satellite of speci f ied a, e, I , 
the minimum adjustment of the semi -major axis is made and the 
nodal longitude is calculated so as to assure optimum cover- 
age of the point mass and sea level test area. For the stan- 
dard case, one close satellite and three distant satellites, 
all geosynchronous, were assumed. Their elements 


a/R 

e 

I 


1.413 

0.01 

90° 


6.61 

0.00 

0° 

212.8° 

6.61 

0.00 

0° 

248.8° 

6.61 

0.00 

0° 

284.8° 


The altitude of the close satellite is thus 263 km. 

3. T racking 

Parameters affecting the entire program are the du- 
ration of tracking and the relative weights of different tracking 
types. The duration used should be the minimum to get good cover- 
age of all elemental blocks in the test area. A close satel- 
lite will have a mean motion near 16 cycles/day* If the period 
is not an exact integral fraction of a day, then low latitude 
squares of side length s° will be crossed about s/11 times a 
day. Durations of 0.8 s days were used in the test runs, thus 
giving about seven crossings for each elemental block. For 
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the standard case, this duration was eight days. 

The weights for all tracking types were unity for the 
standard case. 

The stat ion- to- satel 1 i te laser ranging had a rather 
elaborate error model carried over from the lunar ranging error 
analysis, with allowance for components of random error and bias 
dependent on both zenith angle and sun-satellite angle. For 
the earth physics satellite error analysis the bias was omitted 
and the random error sigma in the standard case was taken to be 

a(r) = [l.O 2 + 0 . 1 2 sec 2 z + 0 . 5 2 cos 2 *]^ (3*0 

where z is zenith angle and ^ is sun-satellite angle. The tracking 
intervals assumed in the standard case were 0.5 days for the 
close satellite and 0.889 days for the distant satellite. The 
program determines when the satellite will be closest to the 
zenith for each station within the specified interval, 0.5 days, 
and assumes the range observed at this time. The tracking inter- 
val for geosynchronous satellites was deliberately made a non- 
integer fraction of a day to obtain a greater variety of geo- 
metries. For tracking over N days, the interval was normally 
taken to be N / (N+1 )/8: i,e., enough to get at least nine ranges 

to determine the six orbital constants-of - integrati on . 

The satel 1 i te- to-satel 1 i te range-rate had only a simple 

random error; +0.5 mm/sec was assumed in the standard case. 

The tracking interval At depended on the spacing s between mass 

points in the test area and the distance r from the nearest 
r cm 
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square in the test area 


At 



( 35 ) 


Based on the criteria that the maximum At anywhere corresponds 
to a 45° step along the orbit, and that the maximum At in or 
into the test area corresponds to a step equal to the spacing s, 
values of k and m are 


s k m 

2.5° 0.00716 1.084 

5° 0.00755 1.014 

10° 0.00793 0.944 

15° 0.00810 0.918 


for At in days and r in earth radii. 

' cm 

The satel 1 i te- to-sea altimetry was assumed to have a 
random error sigma of +2 m. It was assumed that on every pass 
over the test area one observation was made on every elemental 
square crossed. 


4. Natural Parameters 


As described above, the long wavelength part of the 
gravity field was assumed to be represented by a set of 16 11*^ 
and 12*^ degree harmonics, and the short wavelength part by a 
set of nine mass points at 10° intervals over the 30° x 30° square 
The departure of mean sea level from an equipotential was assumed 
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to be represented by the mean values for nine corresponding 
10° x 10° squares. 

In the calculation, only one set of orbital constants- 
of- i ntegrat ion for the distant satellite was carried at a time, 
by using the technique of normal matrix partitioning [ Kauia . 

1966, p. 105]. Hence the total number of parameters for the 
standard case was 16 + 2x9 + 2x6 = 46. 

Runs were designated by a numbering system S - Q, where 
S was the mass point & sea level spacing and Q the number of 
a question to be answered by a variation from the standard case. 

If a question had a minor variant, the variation /lumber was fol- 
lowed by a letter. The standard case was designated 10-1. 

The variations, and the corresponding changed input: 

5-2 Closer resolution than 10°. 

Close satellite a/R = 1.0407 

Satel 1 i te- to-satel 1 i te range-rate interval parameters 
k = 0.00755, m = 1.014 
2.5-2 Close satellite a/R = 1.0404. 

Internal parameters: 

k - 0.00716, m = 1.084 
10-3 Changed tracking station array. 

A. eight stations, two Pacific stations replaced 
by near polar stations: 


20° N, 

204° E 

by 

64° N, 

2 1 2° E : 

Alaska 

26° S, 

216° E 

by 

78° S, 

195° E: 

Antarctica 


B. four stations, three most distant from the 
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test area omitted, leaving: 

38° N, 23° E 
35° N, 202° E 
33° S, 18° E 
33° S, 290° E 

C. 0 stations tracking close satellites -- but 
all eight tracking distant satellites. The corresponding weight 
is set zero. 

10-4 Omitting the satel 1 i te- to-satel 1 i te range-rate. 
The corresponding weight is set zero. 

10-5 Omitting the altimetry. The corresponding weight 
is set zero. 

10-6 Tracking station positions uncertain. 

A. tracking station coordinates are added to 
the parameters set, and assumed to have a priori sigmas of +30 m 

B. the same as 6A, plus a second close satel- 
lite of inclination 45° 

10-7 Only one geosynchronous satellite. The number 
of distant satellites is reduced to one at a longitude over 
the test area. 

10-8 A second close satellite, inclination 45°. 

10-9 Tracking distant satellites at exactly daily 
intervals, instead of 8/9 per day. 
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10-10 Improvement in accuracy of sate) 1 i te-to-satel 1 i te 
range-rate and altimetry by an order of magnitude. 

10-11 Oistant satellites are four in polar orbit in- 
stead of three in geosynchronous equatorial orbit 

10-12 Altitude of close satellites is varied: 

A. lowered by 60 km 

B. rai sed by 60 km 

10-13 Test area at high latitude: 80° +50°. 

10-14 Test area at high latitude and four distant polar 
satellites replace three equatorial geosynchronous satellites. 

10-15 A. One distant satellite, in an orbit such that 
it drifts 360° in eight days: a = 62 1 R . 

B. Three distant satellites drifting 360° in 

eight days. 


IV ALGORITHM 

Before giving the results of computations carried out with 
the standard case and variations therefrom, a description of the 
program procedure may help to understand the nature of the analysis. 

The program has five main components: a main program for 

read-in and conversion of data, matrix manipulation common to 
all derivation types, and write-out of the results; plus four 
principal subroutines, one for each observation type: station- 

to-distant satellite range (GRD); satel 1 i te-to-satel 1 i te range- 
rate (DDC); stat ion- to-c lose satellite range (GRC); and close 
satel 1 i te-to-sea altimetry (CAG). 

The main program reads in all input data, and converts it 
to "planetary" units: i.e., such that the earth's radius and 
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mass and the gravitational constant are unity. All angles are 
converted to radians, and rotation matrices used repeatedly 
are calculated. About the only preliminary adjustments of in- 
put data are to correct the nodal longitude and semi -major axis 
of the close satellites so that the test area is sampled as 

much as possible. The motion of perigee id and the perturbation 

* 

of the mean motion aM by the oblateness J 2 are calculated in 
the usual manner f Kaula . 1966, p. 39]. Luni-solar contributions 
are also included. Then defining 

k = integer part of [ (n + AM + <d )/0 + £] (36) 

the optimized semi -major axis is 

a = u 1/3 [(k - s/2rr)0 - AM- «)]‘ 2/3 (37) 

where 0 is the sidereal rotation rate of the earth and s is 
the size of mean sea level elements in radians. To start the 
satellite off over the test area, the argument of perigee id 
is calculated by 

sin id = sin cp s ysin I, -tt/ 2 < u> < tt/2 (38) 

and the longitude of the node is calculated by 



26 


cos a = cos uj/cos cp 

(39) 

A = + s/2 + 0 q - a 


0 

Figure 2 

where is the longitude of the western border of the test 
area, cp g is the latitude of the south border, and 0 Q is the 
Greenwich sidereal time at the starting epoch. The quadrant 
of a is the same as that of oa. 

The heart of the program is a nest of two "DO" loops, the 
outer over all distant satellites and the inner over all close 
satellites. In the outer, the subroutine for station-to-di stant 
satellite range (GRD) is called first. Then the inner loop 
over close satellites is executed, and for each close satel- 
lite the satel 1 i te- to-satel 1 i te range-rate (DDC), station-to- 
satellite range (GRC) and sate 1 1 i ter to-surface altimetry (CAG) 
subroutines are called. At the end of the outer loop the normal 
equation block corresponding to the distant satellite is reversed 
by the partitioned matrix technique [ Kaula . 1966, p. 105]. 

In the case of the tracking type subroutines (GRD, DDC, 

GRC, CAG) the entire orbit histories are run through for the 
prescribed duration, determining when observations are possible, 
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forming the observation equations at these times, and incre- 
menting the normal equations. 

For station- to-satel 1 i te ranging, the possibility of ranging 
depends not only on geometrical conditions, but also on clear 
weather, as determined by selection of a random number between 
0.0 and 1.0. If this number exceeds the percent clear days, 
then that station does not observe that day. Determination 
of the combination of Greenwich sidereal time 9 and eccentric 
anomaly E for observability is by the iterative technique in 
Kau 1 a [1966, pp. 87-881. The formation of partial derivatives 
with respect to the orbit constants-of- integration and spheri- 
cal harmonic coefficients of the gravitational field is in ac- 
cord with Kau la [1966, pp. 63 & 67-71]. Partial derivatives 
with respect to the mass points are neglected, since the ranging 
will be much less sensitive thereto than are the satel 1 i te- to- 
satel 1 i te range-rate and the altimetry. 

For satel 1 i te- to-satel 1 i te range-rate, after initial cal- 
culations (such as rectangular coordinates of mass points) a 
starting time At is calculated from (35) using the minimum pos- 
sible r cm » a-R. Then intervisibility of the distant and close 
satellites is calculated as the condition that the angle iji sub- 
tended at the earth's center between the two satellites is less 
than the angle * corresponding to a grazing 1 ine-of-s i ght (see 
Figure 3) 
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COS $ > COS t(l, 



cos $ = 


2 2 2 
r c + r d - r l 


cos 


_ f c n d-( r l tr i) i 


2r c r d 


r t > 


[tf -" 2 )*^-" 2 )*] 2 




(40) 


If the satellites satisfy condition (40), the acceleration 
and its derivatives are calculated by (7)- (16). Partial de- 
rivatives with respect to constants-of- integration of both satel- 
lite orbits, spherical harmonic coefficients of the gravitational 
field, and mass points within a maximum range (usually about 

0.4R) are calculated. The minimum r is determined for the 

cm 

purpose of computing the time increment At to the next step’ 
by (35) and the cycle repeated until the end time is reached. 

For satel 1 i te- to-sea altimetry, the mean anomaly and Greenwich 
sidereal time when the satellite will be over the test area 
must be calculated. Let and a ^ be calculated in the same 
manner as a in (38)- (39), but using the latitudes cp^ of each 
of the latitudinal bands of elements in the test area. The 
times of initial crossings, both southward and northward, of 
the central latitude cp_ of the test area are calculated by (refer 
to Figure 2) 
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s 1 n “ Nc = s. 


n ^cj stn -tt/ 2 < u^ < tt/2 


u sc = n • U N 


f = U., - tl) 

N,s N,sc o 


,.„ c c - l/i (l-e 2 ; \ 

cos E N,s " eV “t +e cos f N J 


M N,s * E N,s ' e sin E N,s 


‘N.s “ £ o + K,s- M o)/ A 


(41) 


The node s and Greenwich sidereal times 0 N s are then cal- 
culated, and the longitudes when crossing the central latitude 


X N,sc = a N, sc + n N,s " 9 N,s 


(42) 


If the longitude is within a band making it possible for any 
of the test area elements to be crossed, then the mean anomalies 
for each latitude band are calculated and the same process 
(41)- (42) repeated to determine the exact elements crossed. 

The partial derivatives of altitude with respect to mean sea 

level are calculated by (33). After taking the other partial 

* 

derivatives (29)- (32) and incrementing the normal equations. 
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• e 

the time is advanced by 2tt/(M + <u) and the process repeated un- 
til the whole duration is completed. 

After the loops over all close and distant satellites are 
completed, the main program adds the inverse squares of the 
a priori sigmas to the normal matrix main diagonal; inverts 
the normal matrix; calculates the resulting standard deviations 
and correlation matrix; and prints out these results. 

V RESULTS 


The complete output for the standard case, 10- IA, is given 
on the next five pages. The output sigmas only for all cases 
are given in Table I . The test area was assumed to be between 
longitudes 310°E and 340° E and latitudes 0° and 30° N. The odd 
odd longitudes of the distant satellites given for the stan- 
dard case, 212.8°, 248.8°, and 284.8°, are the consequence of 
an erroneous omission of degree- to-radian conversion which was 
not realized until all the cases had been run. A re-run was 
made with the corrected program, which is designated as 10- IB 
in Table I. The locations of the geosynchronous satellites 
with respect to the test area for 10- 1 A and 10- IB are shown 
in Figure 4. The IB set given one distant satellite directly 
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TABLE I 

Standard Deviations for Parameters of the Gravitational Field Resulting from Error Analyses 


See pp. 7 

- 10 

for def in i t ion 

of parameters. 

and 

pp. 22 - 

24 for 

definition of 

runs . 













Variable 

Un i t 

10- IA 

10-IB 

5-2 

2.5-2 

10- 3A 

10-3B 

10-3C 

10-4 

10-5 

10- 6A 

10- 6B 

10-7 

10-8 

10-9 

10-10 

10- 1 1 

10- 12A 

10-12B 

10-13 

10-14 

10-15A 

10-15B 

c n.o 

10' 9 

.002 

.002 

.002 

.002 

.002 

.003 

.004 

.045 

.002 

.004 

.001 

.002 

.001 

.002 

.0004 

.002 

.001 

.002 

.002 

.002 

.002 

.002 

*'11.3 

<1 

,1*8 

.46 

.71 

1.06 

.48 

.49 

.44 

7.61 

.48 

.49 

.25 

.81 

.25 

.49 

.05 

.29 

.42 

.53 

.48 

.28 

.78 

.45 

S 1 1 , 3 

II 

.28 

.27 

.19 

. 14 

.23 

.34 

.44 

.98 

.28 

.39 

.15 

.37 

.12 

.28 

.05 

.22 

.22 

.24 

.27 

.22 

.42 

.27 

c n ,6 

II 

.63 

.61 

.56 

.60 

.63 

.63 

.60 

10.36 

.63 

.64 

.29 

1.14 

.29 

.63 

.06 

.41 

.55 

.69 

.42 

.25 

1.10 

.65 

S 11,6 

II 

.29 

.27 

. 14 

.06 

.28 

.44 

.59 

.81 

.28 

.51 

.16 

.36 

.12 

.28 

.06 

.22 

.16 

.22 

.20 

.19 

.39 

.27 

C 1 1 ,9 

II 

.39 

.36 

.58 

.91 

.39 

.39 

.35 

7.44 

.39 

.39 

.31 

.63 

.31 

.39 

.04 

.23 

.33 

.42 

.37 

.24 

.64 

.37 

s 1 1 ,9 

II 

‘ .15 

.15 

.13 

.09 

. 13 

.26 

.37 

.68 

.15 

.30 

.15 

.18 

.11* 

.15 

.04 

.13 

. 14 

.17 

.23 

.18 

.20 

.15 

C 1 2 . 0 

II 

.06 

.05 

.02 

.01 

.04 

.09 

.52 

.92 

.06 

.19 

.04 

.06 

.33 

.06 

.03 

.05 

.04 

.07 

.07 

.08 

.07 

.06 

C 12,3 

ll 

.92 

.94 

1.04 

1.03 

.94 

.98 

.92 

8.68 

.92 

.98 

.70 

1.65 

.67 

.94 

.10 

.70 

.76 

.97 

.77 

.51 

1.69 

.93 

S 1 2 , 3 

II 

.62 

.63 

.46 

.20 

.49 

.82 

.93 

1.32 

.62 

.89 

.42 

.79 

.37 

.62 

.10 

.51 

.36 

.70 

.58 

.47 

.83 

.60 

C 12 , 6 

" 

.67 

.63 

.72 

.84 

.67 

.69 

.63 

17.24 

.67 

.69 

.48 

1.19 

.47 

.67 

.07 

.37 

.54 

.69 

.51 

.29 

1.10 

.65 

S 1 2 , 6 

ll 

.44 

.39 

.23 

. 10 

.30 

.59 

.65 

1.08 

.43 

.67 

.25 

.48 

.20 

.43 

.07 

.29 

.31 

.53 

.37 

.25 

.46 

.41 

c 12,9 

ll 

.59 

.51 

.68 

.89 

.59 

.59 

.54 

15.70 

.59 

.59 

.25 

.96 

.25 

.59 

.06 

.30 

.48 

.59 

.58 

.32 

.88 

.52 

S 12,9 

II 

.22 

.20 

.14 

.10 

.20 

.35 

.54 

.42 

.22 

.51 

.15 

.26 

.10 

.22 

.06 

.17 

.21 

.30 

.28 

.24 

.24 

.21 

C 1 2 , 12 

ll 

. 10 

. 10 

.16 

.26 

. 10 

.10 

. 10 

5.94 

.10 

.11 

.10 

.18 

.10 

.10 

.01 

.06 

.09 

.11 

.11 

.06 

.17 

.10 

S 12, 12 

II 

.05 

.05 

.03 

.03 

.04 

.06 

. 10 

. 14 

.05 

.08 

.05 

.06 

.04 

.05 

.01 

.04 

.03 

.04 

.05 

.04 

.06 

.05 

m-NW 

10~ 9 M 

E -55 

.71 

.58 

.64 

.54 

.55 

.52 

25.38 

.55 

.55 

.39 

.81 

.39 

.55 

.06 

.64 

.41 

.58 



1.43 

.56 





















.32 

.38 



m-NC 

II 

.98 

.89 

2.09 

1.32 

.98 

.99 

.86 

34.48 

.99 

.99 

.45 

1.35 

.45 

.99 

. 10 

.88 

.54 

1.01 



2.89 

.96 



.41 


.82 
















.39 

.25 



m-NE 


.35 

1.23 

.41 

.42 

.46 

35.63 

.42 

.42 

.32 

.76 

.32 

.41 

.04 

.40 

.59 

1.04 



.62 

.33 

m-CW 

II 

.41 

.58 

.52 

.64 

.41 

.42 

.40 

45.50 

.41 

.42 

.32 

.77 

.32 

.42 

.04 

.49 

.27 

.57 

.48 

.67 

1.14 

.42 

m-CC 

•1 

- .74 

.65 

2.35 

1.35 

.74 

.74 

.78 

47.31 

.74 

.74 

.40 

1.29 

.40 

.74 

.07 

.74 

.51 

1.00 

.42 

.48 

1.11 

.77 

m-CE 

II 

.37 

.28 

.89 

1.25 

.37 

.38 

.44 

45.69 

.37 

.38 

.30 

.55 

.30 

.38 

.04 

.47 

.67 

.98 

.44 

.42 

.49 

.29 

m-SW 

II 

.47 

.56 

.48 

.62 

.47 

.47 

.34 

34.60 

.47 

.47 

.36 

2.15 

.36 

.47 

.05 

.44 

.22 

.64 

.84 

1.16 

4.29 

.48 

m-SC 

II 

.73 

.70 

1.86 

1.27 

.73 

.74 

.74 

35.44 

.73 

.74 

.47 

1.42 

.47 

.75 

.08 

.68 

.59 

1.07 

.79 

.80 

1.10 

.76 

m-SE 

II 

.40 

.29 

.81 

1.17 

.40 

.40 

.42 

28.09 

.40 

.40 

.32 

.72 

.32 

.41 

.04 

.71 

.94 

1.76 

.78 

.65 

.46 

.31 
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over the test area, and a second which can see the close satel- 
lite when it is over the eastern edge. The IA set gives two 
distant satellites off to the side, subtending an angle of about 
38° at the test area, and seeing the close satellite anywhere 
over the area. The consequence appears to be that IB is rela- 
tively better for the eastern edge and I A is better for the 
western edge. Since none of the sigmas differ by more than 
25%, it was decided that it was not worthwhile re-running all 
cases with the corrected longitude (except 10-7, a single satel- 
lite over the test area). 

The results from runs 5-2 and 2.5-2 indicate the resolution 

limit of the +0.5 mm/sec sate 1 1 i te- to- sate 1 1 i te range-rate with 

the 260 km close satellite altitude is somewhere around 3° or 

4°: say 350 - 400 km. The 5° output sigmas for the point masses, 

0.48 to 2.35 x 10 , are wel 1 under the a priori , 6. Ox 10 (see 

p. 9). The 2.5° output has six sigmas which are only slightly 

-9 

reduced from the a priori; 1.5 x 10 , even though the calcu- 

lation was carried twice as long as for 5-2 and four times as 
long as for 10- 1 A. This result is somewhat affected by the 
imperfection of the algorithm for determining the semi -major 
axis, described on pp. 25- 26; for case 2.5-2, the area was 
not uniformly covered. 

The results for 10-3A suggest that tracking stations near 
the pole are helpful for determining the low order (small m) 
harmonics. The results for 10-3C, 10-4, and 10-5 suggest that 
both the +0.5 m ranging from tracking stations to the close 
satellite and the +2 m altimetry have little weight in deter- 
mining the gravity field variations compared to the +0.5 mm/sec 
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satel I i te- to-satel 1 i te range-rate. The large increase q) 

from 10- 1 A to 10-3C suggests that the principal value of the 
ranging from ground stations is to determine the orientation 
of the orbit. This is corroborated by the increases in a^Q c j 
and ct(uu c ) between the two runs, which were by factors of 3 and 
respectively. 

The results for 10-6A indicate that tracking station errors 
have some effect on determining the broad variations of the 
gravity field, but not the shorter wavelength variations here 
represented by mass points. This result is somewhat affected 
by the interval between satel 1 i te- to-satel 1 i te range rates being 
as high as 45° far from the test area. 

The poorer results from 10-7, on which a single geostationary 
satellite over the area was used, indicate that viewing from 
more than one direction is helpful in resolving the shorter 
wavelength variations in the gravity field. The results from 
1 0- 1 5 A , the single drifting distant satellite, are disappointing 
as a means of obtaining this variety of directions; possibly 
the eight day duration of the test is too short -- certainly 
the geostat ionary satellite locations with respect to the test 
area in 10- IA are better than the average expected. 

The results from 10-11 indicate that tracking from the 
equatorial geosynchronous distant satellites can resolve variations 
in the gravity field at high latitudes as well as at low. The 
results from 10-13 and 10-14 further indicate that the four 
distant satellites in polar orbits are not significantly better 
or worse than the three equatorial geosynchronous satellites. 
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The second close satellite at a lower inclination, tested 
by runs 10-6B and 10-8, results in the expected improvement: 
i.e., the sigmas are reduced by a factor of about 1//?. Some 
correlations between gravity coefficients are reduced from about 
0.6 to 0.3 or 0.4 (e.g., ^ and Cj 2 0 and C 12 9)* 

The distant satellite tracking interval, tested by run 
10-9, has a negligible effect on the gravity field determination, 
but it did have some effect on determination of distant satel- 
lite orbit elements, the largest increase being in a(e), by 
a factor of 2. It would also have some effect on the determi- 
nation of station location. 

Significant determination of ri 0 j , the mean sea level de- 
parture from the geoid, is obtained only when much better alti- 
metry is assumed, run 10-10. From the +10 cm altimetry over 
eight days sigmas of +4 to 6 cm were obtained. 

Variation of the close satellite orbit altitude by +23% 
in runs 10-12A and 1 0- 1 2B produced appreciable changes in the 
sigmas for the point masses, on the order of +50%. Some of 
the changes for the lower altitude satellite were not in the 
expected direction, apparently because the orbital node and 
semi-major axis could not be optimized. But it seems clear 
the message is to keep the close satellite as low as practicable. 

Comparison with Other Results 

Error analyses of sate 1 1 i te- to- sate 1 1 i te tracking have 
also been done by Martin et al . [1972] and Schwarz [1972], 
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That by Hart in et a 1 . [1972] assumed range, rather than range- 
rate, between the satellites, and analyzed only the determi- 
nation of low degree spherical harmonic coefficients. Hence 
a meaningful comparison is difficult. 

The analysis by Schwarz [1972] assumed range-rate tracking. 
Values of the surface density coating were assumed for 5° x 5° 
and 2° x 2° rectangles and solutions were made from simulations 
of orbital passes and associated tracking. Both distant-close 
and close-close satellite pairs were tested. It was concluded 
that the technique could resolve elements about the same size 
as the satellite altitude, the definition of resolution being 
that the correlation coefficients between adjacent blocks are 
on the order of 0.8. In the present study, the correlation 
coefficients between adjacent blocks were about 0.5; however, 
since the output sigmas approached the a priori sigmas, this 
value is probably lower than would have been attained without 
a priori sigmas. To the extent that the two analyses can be 
compared, the agreement seems satisfactory. 

VI CONCLUSIONS 
Principal Recommendations 

The most useful tool for determination of the variations 
in the gravity field appears to be definitely the satellite- 
to-satellite tracking (unless and until some sort of gradio- 
meter proves feasible); +0.5 nm accuracy appears capable of 
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resolving features a little greater in extent than the satel- 
lite altitude. The orbit should be kept as low as practicable; 
some sort of orbit maintenance system should be considered. 
However, since satel 1 i te-to-satel 1 i te tracking requires determi- 
nation of a reference orbit, the system should be operated inter- 
mittently, not frequently. The variations in accelerations 
due to drag should be small compared to the gravitational varia- 
tions of comparable wavelength at any altitude at which an orbit 
can be maintained. There appears to be a wide variety of allow- 
able geometries of distant satellites, provided that the close 
satellite is tracked from at least two distinctly different 
directions; whether this requirement can be satisfied by tracking 
at different times is not yet clear. 

Satel 1 i te-to-satel 1 i te tracking would definitely allow 
economies of tracking from ground stations. How much is not 
clear from the present study, because the intervals of satel 1 i te- 
to-satel lite tracking were deliberately varied with distance 
from an assumed test area. Tracking of distant satellites from 
ground stations remains necessary, of course. 

The altimeter seems inherently less capable of resolving 
details in the gravity field, because it measures an integral 
of what is measured by the satel 1 i te-to-satel 1 i te range-rate. 
However, since the present error of the geoid is about +10 m rms, 
a +2 m system would manifestly gain considerable information. 

To resolve differences between geoid and sea level accuracies 
on the order of +20 cm are needed. 


■y 
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Limitations and Possible Future Work 

In modeling the satellite and tracking systems, the pre- 
sent program needs to be extended to include range-rate between 
two close satellites, and VLBI to distant satellites. Instru- 
mental biases need to be examined with a view to more elabo- 
rate modelling, such as has been applied in error analyses of 
lunar ranging [ Kaula . 1972], More investigation of alternative 
geometries of orbits than was done in the present study should 
be undertaken. 

In the modeling of natural parameters, there still re- 
main several deficiencies. Time varying effects on the sea 
level, both regular (tides) and stochastic (storms) should be 
included. Because of the broad spectrum entailed, it would 
take some thought how to do this economically. The short wave- 
length variations in the gravity field should be represented 
by surface coatings for area elements, rather than by point 
masses. The perturbations of the radial coordinates of the 
close satellite by short wavelength variations should be in- 
cluded; how to do this with an orbital representation depen- 
dent on a mean intermediary is not yet clear. 
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